Statistical Mechanics of a Discrete Schrodinger Equation with Saturable Nonlinearity 
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We study the statistical mechanics of the one-dimensional discrete nonlinear Schrodinger (DNLS) 
equation with saturable nonlinearity. Our study represents an extension of earlier work [Phys. Rev. 
Lett. 84, 3740 (2000)] regarding the statistical mechanics of the one-dimensional DNLS equation 
with a cubic nonlinearity. As in this earlier study we identify the spontaneous creation of localized 
excitations with a discontinuity in the partition function. The fact that this phenomenon is retained 
in the saturable DNLS is non-trivial, since in contrast to the cubic DNLS whose nonlinear character 
is enhanced as the excitation amplitude increases, the saturable DNLS in fact becomes increasingly 
linear as the excitation amplitude increases. We explore the nonlinear dynamics of this phenomenon 
by direct numerical simulations. 

PACS numbers: 05.45.-a, 63.20.Pw, 63.70.+h 

During the past two decades much research has been 
devoted to nonlinear mechanisms for the storage and 
transport of localized coherent packages of energy and 
charge in spatially periodic systems, which may be de- 
scribed by discrete lattice models in one, two, or three 
spatial dimensions [TH3]. Generally such systems may 
support intrinsic localized modes, or discrete breathers, 
which are time-periodic, spatially localized solutions to 
the dynamical lattice equations [4 . 

A much older interest [5] in such nonlinear systems 
arises from their ability to thermalize through their in- 
trinsic dynamics leading to equipartition of excitation en- 
ergy among the linear modes. Localization and equiparti- 
tion appears somewhat contradictory and we have earlier 
explored their interrelation within the framework of the 
Discrete Nonlinear Schrodinger (DNLS) equation with a 
cubic nonlinearity [6} Cj . We found that in the case of the 
cubic DNLS these phenomena occur in distinct param- 
eter domains. Thermalization in the Gibbsian sense oc- 
curs in the low energy part of the parameter space while 
the nonlinear dynamics for higher energies lead to the 
spontaneous creation of localized excitations prohibiting 
thermalization. It has been argued (8] that the statistical 
process behind this localization leads towards an entropy 
maximization. The entropy maximization is driven by 
optimally distributing the effects of fluctuations (which 
carry the main part of the system entropy) among dif- 
ferent conserved quantities. Localization is therefore ab- 
sent (Gibbsian regime) if insufficient energy is supplied 
by the initial conditions. If a surplus of energy is pro- 
vided by the initial conditions, it is allotted to the high 
amplitude structures (non-Gibbsian regime) that absorb 
large amounts of energy while using few lattice positions. 
The specific separation of these two regimes was derived 
analytically [B]. 

Extending this prior work we will here detail a sim- 
ilar study for the DNLS with a saturable nonlinearity, 
which physically describes, e.g., discrete spatial solitons 
in a tight-binding approximation of one-dimensional op- 



tical waveguide arrays made from photorefractive crys- 
tals [9]- Another example where this equation arises is 
that of Bose-Einstein condensates in optical lattices [TU] . 
Mathematically the saturable DNLS reduces to the cubic 
DNLS in the low amplitude (energy) limit, while it be- 
comes essentially linear in the large amplitude limit. The 
saturable DNLS is known to support localized excitations 
[TTHHj, but it is unclear how the interplay of thermal- 
ization and spontaneous localization occurs in this case 
where strong localization essentially leads to linear dy- 
namics. 

To explore these issues we are concerned with the 
DNLS equation with a saturable nonlinearity 
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or equivalently through the simple gauge transformation 
<t>n = i>n exp(ii/t) 
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where ip n (4> n ) is a complex valued "wave function" at 
site n and v is a real parameter that controls the non- 
linearity of the system. It is important to note that for 
small excitation amplitudes \vb n \ <c 1 (\<t>n\ "C 1) both 
equations reduce to the cubic DNLS. 

Equation (f2j) represents a Hamiltonian system for the 
canonically conjugated pairs iip n and -0* with the Hamil- 
tonian 
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In addition 



so that Eq. &2i is given by iip n — 

to the Hamiltonian, H, the quantity A = J2n=i IVvJ 2 is 
also conserved by the dynamics of Eq. pi) and serves as 
the norm of the system. 



Our objective is to calculate the grand-canonical par- 
tition function Z of this system and to accomplish this 
we first apply the canonical transformation \ip n \ 2 — 4 n , 
tpn = \/A n e l6n , which brings the Hamiltonian % of Eq. 
([3]) into the form 
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and the norm becomes A = J2 n =x A n- 

In terms of these variables the partition function Z is 
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where the parameter /x is introduced analogous to a chem- 
ical potential to ensure the conservation of A. From this 
partition function all thermodynamic quantities of the 
system can be calculated. We will in particular be inter- 
ested in the averaged norm density a — (A)/N and the 
averaged energy density h — {%) /N , which are given by 
the following expressions 

1 dlxiZ , 1 dlxiZ 

a = — ^ttt — ^ — ; h + fia = — — aQ . (6) 
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In Eq. ([5]) the integration over the phase variables 9 n 
can be preformed straightforwardly reducing the parti- 
tion function to 
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where Iq denotes the modified Bessel function. Further 
treatment of the partition function has to be carried out 
numerically using for example the transfer integral ap- 
proach [5J. However, in the two limits p — > 00 and p — > 
we can gain further analytical insight. Noticing that the 
Hamiltonian is bounded from below one can observe that 
its minimum (corresponding to P — > 00) is realized by 
the plane wave ip n = ^/a exp imr , whose energy density 
is ft. = —2a — v ln(l + a). When (3 —>■ while /3/i remains 
finite the modified Bessel function in the expression for 
Z can reasonably be approximated by unity, which leads 
to 

N »oo 
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Making the substitution 1 + A n = x n brings the partition 
function Z into 

Z = {2ir) N e N ^E(pv, p^) N , (9) 

where the function E(x, y) is given by 
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Keeping Pfi = 7 finite as we take the limits P — > and 
/1 4 00, we can, using Eq. ([6]) and the properties for 
E(x, y) listed in the appendix, obtain expressions for a 
and h 
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and 
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so that 
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We can now compare with the result of the cubic DNLS 
by considering small a. Using the last equality of Eq. 
(A.4| in the limit of small a, Eq. (14 1 becomes 
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In Fig. [T] we illustrate, for v = 1, these results in 
the (a, h)-p\ane. The thick lines represent /3 = [Eq. 
( Ill )] and /3 = 00 [/i = -2a - t/ln(l + a)] that bound 
the regime in which Gibbsian thermodynamics applies. 
The region below the minimum energy line p = 00, is 
of course inaccessible and is therefore shaded. However, 
the region above the /3 = line is accessible but can- 
not be captured by Gibbsian thermodynamics. In the 
corresponding region for the cubic DNLS we found [BJ 
that the nonlinear dynamics lead to the spontaneous ap- 
pearance of strongly localized discrete breather excita- 
tions. The dots overlaid on the P = line represent 
values obtained numerically by applying the transfer in- 
tegral technique (see details in Ref. |BJ) to Eq. ^ using 
P = 10~ 8 . The numerical solution does not ignore the 
modified Bessel function in Eq. (fsj and it therefore di- 
rectly verifies the validity of the expression in Eq. (fl4J) . 
Finally, the dashed line illustrates the small amplitude 
approximation in Eq. (15), which corresponds to the 
cubic DNLS [T5]. In order to explore the differences in 
the dynamics in the two regions separated by the /3 = 
line we perform a number of direct numerical simula- 
tions of Eq.S, with v = 1. These numerical simulations 
represent a micro-canonical ensemble since they are per- 
formed to conserve both norm and energy. We use initial 
conditions of the form ip n = y / aexp(i^n), which can be 
realized in both regimes by tuning the wave number 9 
and the norm a. The energy of these initial conditions is 
given by h = 2a cos 9 — ln(l + a). 




FIG. 1. Parameter space (a, h) [for v — 1], where the shaded 
area is inaccessible. The thick lines represent ft = [Eq.( 14 1] 
and f5 = oo [h — —2a — ^ln(l + a)] and thereby bound 
the regime in which Gibbsian thermodynamics applies. The 
points are obtained numerically by applying the transfer in- 
tegral technique to Eq. |8} using /3 = 10 -8 . The dashed line 
represents Eq. ( |15[ ) which corresponds to the cubic DNLS. 
Initial conditions for direct (micro-canonical) numerical cal- 
culations (see Figs. [21 py and B| are marked by a square, a 
circle, an up-pointing triangle, a down-pointing triangle, and 
a star. 



Our choices of initial conditions are marked in Fig. [T] 
and the result of simulations with these initial conditions 
are shown in Figs. [2] and [3] Figure [2] shows the ampli- 
tude distributions p(A) resulting from long-time simula- 
tions. The numerically obtained distributions are shown 
with triangles, circles, and squares while the distribu- 
tion functions obtained by the transfer integral approach 
are shown with thick lines. The thin lines show the 
function exp(— A/a)/a. As the above analytical work 
demonstrates, this function approximates the distribu- 
tion function in the limit of small /?. The validity of this 
approximation is clearly verified by Fig. [2] since the thin 
lines are only discernible for /3 sufficiently large. Figure 
[2] also clearly demonstrates that the distribution arising 
from initial conditions above the f3 = line arc signifi- 
cantly different from the Gibbsian distribution functions. 
While the Gibbsian distribution functions are concave 
the non-Gibbsian distribution functions are convex and 
have significant tails at large values of A. This convex- 
ity is reminiscent of the discussion in Ref. [5] regarding 
"negative temperatures" and is a subject that recently 
has been discussed in more detail |16j in this context. 
These features are even more pronounced in the example 
given in Fig. [3] which originates from initial conditions 
located deep within the non-Gibbsian regime. The long 
tail distributions indicate the presence of a large ampli- 
tude excitations in the dynamics. The abrupt drop in 



FIG. 2. Distributions of A = \ip n \ 2 for four different ener- 
gies at a = 1.5 [for v = 1]. The symbols show data obtained 
by direct numerical solution. Solid lines show the distribu- 
tion function obtained by solution (using the transfer integral 
technique) of Eq. (ISj) , while the thin lines show the function 
exp(— A/a)/a. The curves have been shifted vertically to fa- 
cilitate visualization. All data is for a — 1.5 while the energies 
are h — —2.13 (squares), h = —0.97 (circles), h — —0.76 (up- 
pointing triangles), and h = 0.70 (down-pointing triangles). 
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FIG. 3. Distributions of A = |V>n| 2 for h = 2.09 at a = 1.5 
[for v — 1]. The distribution is obtained by direct numerical 
solution of Eq. pi). 



the distribution function at very small excitation ampli- 
tudes clearly shows that the large amplitude excitations 
exist on an extensive background of small amplitude ex- 
citations. This scenario is quite similar to Rumpf's [8 
description of the dynamics of the cubic DNLS in this 
regime, which in turn suggests that also in this case does 
localized high-amplitude excitations absorb a surplus of 



energy when they emerge as a consequence of the pro- 
duction of entropy in the small uctuations. 

In Fig. [4] we directly illustrate these excitations with 
a spatio-temporal snapshot of the dynamics. A few very 
large amplitude excitations occur early in the dynam- 
ics. These excitations are much larger in amplitude than 
what is typically observed in the cubic DNLS and al- 
though they are very discrete, only spanning 1-3 sites 
they are highly mobile [14] . This mobility allows the ex- 
citations to collide and these collisions appear to lead to 
further amplitude enhancement. In spite of this high mo- 
bility the probability distribution function still reaches 
stationarity very slowly and require integration times 
above 10 6 time units. The process remains slow because 
the collisions become quite rare as the excitation ampli- 
tudes increase. 
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FIG. 4. (Color online) Spatio-temporal snapshot of the dy- 
namics in the strongly non-Gibbssian regime for h = 2.09 at 
a = 1.5 [for v = 1]. 



In summary, we have studied a nonlinear Schrddinger 
equation with saturable nonlinearity, and derived the sep- 
aration between Gibbsian and non-Gibbsian regime in 
the (a, h) phase space. Through analytical calculations 
supported by direct numerical simulations wc have, as in 
the case of the cubic DNLS, been able to link the non- 
Gibbsian regime to the appearance of localized modes. 
The fact that this phenomenon is retained in the sat- 
urable DNLS is significant since in contrast to the cubic 
DNLS whose nonlinear character is enhanced as the ex- 
citation amplitude increases the saturable DNLS in fact 



becomes increasingly linear as the excitation amplitude 
increases. 
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Appendix: Properties of E(x, y) 

Some specific properties of the function E(x, y) and its 
derivatives: 



E(x,y) 
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For the n'th derivative of E(x,y) with respect to x we 
find 
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and for the n'th derivative of E(x,y) with respect to y 
we find 
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Using these expressions we find the following 
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The last expression follows for example from Eq. (5.1.51) 
of Ref. [T7] where the exponential integral is defined as 

E x {y) = J™ ^-dt (see Eq. 5.1.1 of Ref. [17). 
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